Introduction

Cruciferous vegetables, including broccoli, Brussels sprouts, collard greens, and arugula, ubiquitously contain a class of compounds known as glucosinolates (GLS) [1]. A single type of cruciferous vegetable will contain a wide variety of GLS, however, they typically contain one to three GLS in the greatest abundance [2]. Following damage to the cell wall of the vegetable, the enzyme myrosinase converts GLS to bioactive isothiocy-anates (ITC), or epithiospecifier protein can facilitate conversion to biologically inert ni-triles (NIT) [1,2]. Broccoli sprouts are known to be rich in the aliphatic GLS glucoraphanin and also contain glucoiberin which yield the ITCs sulforaphane (SFN) and iberin (IBN), respectively [2]. Epidemiological studies have shown an inverse relationship between the incidence of lung, breast, prostate, colon, and bladder cancer and cruciferous vegetable consumption, offering an appealing, cost-effective, non-pharmacological approach to cancer prevention through dietary intervention [3–6]. Furthermore, cell culture and animal models have shown that aliphatic ITCs, such as SFN, have potent anticancer properties and can suppress and prevent the formation of cancer [7–16] . In contrast to these findings, human clinical trials examining the efficacy of cruciferous vegetables interventions on cancer prevention targets have shown high levels of inter-individual variation in both the absorption and excretion of ITCs, however, the source of this variation is unknown (reviewed in [17]). Early work has focused on glutathione-S-transferase polymorphisms, as a major factor modulating effects of ITC metabolism and variable effects on cancer, but studies are inconsistent and equivocal [18–30]. Our own group has shown that Nrf2 KO mice, which lack ability to induce glutathione-S-transferases do not manifest differences in SFN metabolite production [31]. As an alternative, it has been proposed that differences in individual gut microbiomes may contribute to the observed variation in ITC levels, potentially through the production of biologically inert NITs from GLS by gut microbes. This is significant as an increase in NIT production from GLS reduces the amount of ITCs produced and thus may lower the efficacy of cruciferous vegetables interventions.

Members of the human commensal gut microbiota have been shown to hydrolyze GLS to both ITCs and NITs in the absence of myrosinase, such as in cooked vegetables [32–35]. Preliminary studies conducted in humans fed glucoraphanin show variations in SFN excretion coincide with variations in gut microbiome composition [36,37]. In addition to this evidence, in vitro studies demonstrated that different microbes preferentially hy-drolyze GLS to either ITCs or NITs [38–41]. Additionally, consumption of cruciferous vegetables can significantly alter the composition of the gut microbiome highlighting potential diet-microbiome effects [36,37,42]. It is hypothesized that the gut microbiome is the underlying driver of variation observed in many human clinical trials, however, evidence in complex and dynamic systems is lacking. More specifically, work on the role of the gut microbiome in GLS metabolism is limited primarily to in vitro studies conducted using GLS-isolates in monocultures of culturable organisms. These studies fail to capture both the impacts of a whole food matrix on the microbial metabolism of GLS, and the effects of broccoli consumption on microbiome composition.

To further investigate the role of the gut microbiome in the production of NIT from GLS, we utilized an ex vivo fecal incubation model using individual human stool samples and broccoli sprouts that underwent in vitro digestion. Despite a growing interest in how the microbiota affects dietary chemoprevention agents, few studies have closely examined the interaction among the gut microbiota and GLS in complex and dynamic systems. The results of this study will assist in understanding determinants of individual response to cruciferous vegetable metabolism and responses to dietary interventions.

Experimental Methods

We completed an in vitro digestion of cruciferous vegetables then incubated them with human fecal stock for 24 hours under anaerobic conditions.

For this experiment we had 3 treatment groups:

  • fecal_stock : The fecal slurry before incubation
  • no_veg : Fecal slurry incubated with media under anaerobic conditions (NC in the paper)
  • broc : 1/2 cup fresh broccoli sprouts in vitro digested then incubated with fecal slurry under anaerobic conditions

Following incubation, untargeted metabolomics was completed in positive and negative ionization modes and 16S rRNA sequencing was competed.

#Coding Tools
library(tidyverse)
library(magrittr)
#Micorbiome Tools
library(phyloseq)
library(vegan)
library(dunn.test)
#library(DESeq2)
#Visualization Tools
library(DT)
library(cowplot)
library(plotly)
library(RColorBrewer)

#Set up color palette:
paltrt <- RColorBrewer::brewer.pal(name = 'Dark2', n = 6)
palsmp <- viridis::viridis(n = 10)
palgrp <- c('#A6611A', '#018571')

#Load in the data
asvtab <- readRDS('./data/asv_tab.RDS')
taxatab <- readRDS('./data/tax_tab.RDS')
metadata <- read.csv('./data/microbiome_metadata.csv')
metabinterest <- read_csv('./data/metab_interest.csv')
rawdata_atwell <- read_csv('./data/atwell_nitriles_rerun.csv')
volumes <- read_csv('./data/atwell_volumes.csv') 

metadata %<>% mutate(treatment = ifelse(treatment == 'no_veg', 'NC', treatment))

#Factorize and set the levels on the metadata
metadata$treatment %<>% factor(levels = c('fecal_stock', 'NC', 'broc', 'brus', 'combo', 'control_digest'))
metadata$fecal_sample %<>% factor(levels = c('T5631','T5632','T6260','T6291','T4669','T1995','T5627','T5717','T5854','T6382')) 

#Set up colors for each treatment
names(paltrt) <- levels(metadata$treatment)
colTrt <- scale_colour_manual(name = "treatment", values = paltrt)

names(palsmp) <- levels(metadata$fecal_sample)
colSmp <- scale_colour_manual(name = "fecal_sample", values = palsmp)

names(palgrp) <- levels(metadata$group)
colGrp <- scale_color_manual(name = "group", values = palgrp)


#Construct the phyloseq object
rownames(metadata) <- metadata$sample
rownames(asvtab) <-metadata$sample
ps_raw <- phyloseq(otu_table(asvtab, taxa_are_rows = FALSE),
               sample_data(metadata),
               tax_table(taxatab))

#Give arbitrary names to the taxa as opposed to keeping as just DNA-sequences which identify them
taxa_names(ps_raw) <- paste0("ASV", seq(ntaxa(ps_raw)))

#Fill in missing genus names:
renames <- rownames(tax_table(ps_raw)[is.na(tax_table(ps_raw)[, 'Genus'])])
taxdf <- tax_table(ps_raw)[renames,]
renamed_genus <- unname(sapply(taxa_names(taxdf), function(x) paste0('f_', taxdf[x, 'Family'], '_', x)))
tax_table(ps_raw)[renames, 'Genus'] <- renamed_genus

#Remove the control digests, these are not relevant to our analysis
ps_raw <- ps_raw %>% subset_samples(treatment != 'control_digest')
#Agglomerate to the genus level
ps_genera <- ps_raw %>% tax_glom(taxrank = "Genus")
#Remove taxa not seen more than 3 times in at least 20% of the samples
ps_counts <- ps_genera %>% filter_taxa(function(x) sum(x > 3) > (0.2*length(x)), TRUE)
#Convert from counts to relative abundance
ps_relab <- ps_counts %>% transform_sample_counts(function(x) x / sum(x) )
#Filter out low abundance (>1e-5) taxa
ps <- ps_relab %>% filter_taxa(function(x) mean(x) > 1e-5, TRUE)

#Number of ASVs of each
summary <- data.frame(
  list('ps_raw' = dim(tax_table(ps_raw))[1],
  'ps_genera' = dim(tax_table(ps_genera))[1],
  'ps_counts' = dim(tax_table(ps_counts))[1],
  'ps_relab' = dim(tax_table(ps_relab))[1],
  'ps' = dim(tax_table(ps))[1]
))
rownames(summary) <- 'Number of ASVs'

knitr::kable(summary, caption = 'ASVs across data preprocessing')
ASVs across data preprocessing
ps_raw ps_genera ps_counts ps_relab ps
Number of ASVs 3557 935 72 72 72 
#Extract relevant sub-data:
ps_stock <- ps %>% subset_samples(treatment == 'fecal_stock')
ps_broc <- ps %>% subset_samples(treatment == 'broc')
ps_bnv <- ps %>% subset_samples(treatment %in% c('NC', 'broc'))
ps_fbnv <- ps %>% subset_samples(treatment %in% c('fecal_stock', 'NC', 'broc'))


datatable(data.frame(sample_data(ps_fbnv)),
           filter = 'top',
           extensions = 'Buttons',
           options = list(pageLength = 12,
                          dom = 'Bfrtip',
                          scrollX = TRUE))

We start by investigating the composition of our pre-incubation fecal slurries

Characteristics of the Pre-Incubation Fecal Slurries

Alpha Diversity

Fecal Stock Composition

We start by investigating the alpha diversity of each of our fecal stocks.

Alpha Diversity Metrics

adiv_stock <- ps_raw %>% 
  subset_samples(treatment == 'fecal_stock') %>%
  estimate_richness(measures = c('Observed','Shannon','Simpson')) %>%
  cbind(., data.frame(sample_data(ps_raw %>% 
                         subset_samples(treatment == 'fecal_stock')))%>%
                         dplyr::select(fecal_sample)) %>%
  remove_rownames() %>%
  column_to_rownames('fecal_sample')

gstock <- adiv_stock %>%
  rownames_to_column('fecal_sample') %>%
  pivot_longer(where(is.numeric), names_to = 'measure')

stk <- ggplot(gstock, aes(x = fecal_sample, y = value, color = fecal_sample)) +
  geom_point() +
  facet_wrap(vars(measure), scales = 'free', ncol = 3) +
  theme_cowplot() +
  theme(axis.text.x = element_text(angle = 90), axis.title.x = element_blank()) 
stk

#apply(adiv_stock , 2, mean)

knitr::kable(adiv_stock, caption = 'Fecal Stock Alpha Diversity Metrics')
Fecal Stock Alpha Diversity Metrics
Observed Shannon Simpson
T4669 578 5.964509 0.9964655
T6382 364 5.281253 0.9913638
T5627 580 5.916242 0.9961635
T5632 441 5.322444 0.9913482
T6260 456 5.535625 0.9934348
T5854 595 6.048985 0.9966105
T5631 696 6.220348 0.9972688
T5717 623 6.174746 0.9973943
T1995 482 5.744007 0.9950483
T6291 678 6.291398 0.9977060

Let’s see if the differences in alpha diversity that we see differ by age or biological sex.

st <- adiv_stock %>%
  rownames_to_column('fecal_sample') %>%
  left_join(.,data.frame(sample_data(ps_stock))) %>%
  select(fecal_sample, Observed, Shannon, Simpson, birth_year, biological_sex) %>%
  mutate(by_bin = ifelse(birth_year <= 1975, 1, 
                         ifelse(birth_year <= 1985, 2, 
                                ifelse(birth_year <= 1995, 3, 4)))) %>%
  pivot_longer(cols = c('Observed', 'Shannon', 'Simpson'), names_to = 'measure') %>%
  group_by(measure) %>%
  nest() %>%
  mutate(kw_b = map(data, function(x) kruskal.test(value ~ as.factor(by_bin), data = x) %>% broom::tidy() %>% pull('p.value'))) %>%
  mutate(kw_s = map(data, function(x) kruskal.test(value ~ biological_sex, data = x) %>% broom::tidy() %>% pull(2))) %>%
  unnest(cols = c('kw_b', 'kw_s'))

stc <- st %>%
  select(-data) %>%
  rename('Biological Sex' = kw_s, 'Age' = kw_b, 'Measure' = measure)

knitr::kable(stc, caption = 'Kruskal Wallis Test Results (p.values)')
Kruskal Wallis Test Results (p.values)
Measure Age Biological Sex
Observed 0.2406619 0.2864220
Shannon 0.2230772 0.2008251
Simpson 0.1856053 0.1355930

We found that there are no significant differences by age nor biological sex. Next let’s investigate beta diversity between our stocks.

Beta Diversity

ordBCstock <- ordinate(ps_stock, method = 'PCoA', distance = 'bray')
sbd <- plot_ordination(ps_stock, ordBCstock, color = 'fecal_sample') +
  geom_point(aes(text = paste0('Fecal Sample: ', fecal_sample, '\n',
                               'Biological Sex: ', biological_sex, '\n',
                               'Year of Bith: ', birth_year)),
             size = 3) +
  theme_cowplot() +
  colSmp +
  ggtitle('Beta Diversity of Fecal Stocks')

ggplotly(sbd, tooltip = 'text', width = 900, height = 600)

Let’s see if the patterns we see in our beta-diversity can be explained by age or biological sex. We will do this via a PERMANOVA test.

PERMANOVA by Age

ordBCstock <- ordinate(ps_stock, method = 'PCoA', distance = 'bray')
sbd <- plot_ordination(ps_stock, ordBCstock, justDF=T) 
distmat <- phyloseq::distance(ps_stock, method = 'bray')
#Extract the metadata
mdata <- data.frame(sample_data(ps_stock)) %>% dplyr::select(fecal_sample, treatment, biological_sex, birth_year)  %>%
  mutate(by_bin = ifelse(birth_year <= 1975, 1, 
                         ifelse(birth_year <= 1985, 2, 
                                ifelse(birth_year <= 1995, 3, 4)))) %>%
  modify_at('by_bin', as.factor)

#Set the seed for reproducibility
set.seed(120)
#Run betadisper to verify the distrubtion of our groups is equal, an underlying assumption of PERMANOVA
bdisp <- betadisper(distmat, mdata$biological_sex)
#Evaluate using a permutation test
permutest(bdisp) 
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.003881 0.0038807 0.1409    999  0.706
## Residuals  8 0.220371 0.0275464
adonis2(distmat~biological_sex, data = mdata)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distmat ~ biological_sex, data = mdata)
##                Df SumOfSqs      R2      F Pr(>F)
## biological_sex  1  0.15232 0.12581 1.1513  0.331
## Residual        8  1.05838 0.87419              
## Total           9  1.21070 1.00000

PERMANOVA by Birth Year

bdisp <- betadisper(distmat, mdata$by_bin)
#Evaluate using a permutation test
permutest(bdisp) 
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq     F N.Perm Pr(>F)
## Groups     3 0.097361 0.032454 1.376    999  0.344
## Residuals  6 0.141510 0.023585
adonis2(distmat~by_bin, data = mdata)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distmat ~ by_bin, data = mdata)
##          Df SumOfSqs     R2      F Pr(>F)
## by_bin    3   0.3781 0.3123 0.9082  0.537
## Residual  6   0.8326 0.6877              
## Total     9   1.2107 1.0000

Neither age (binned in 10-year intervals) nor biological sex were found to be significant as determined by a permutation analysis of variance (PERMANOVA) test, indicating that these factors were not driving differences in microbial composition between fecal stocks.

Incubation with Broccoli Sprout Digest Does Not Alter Gut Microbial Composition

Next let’s investigate the impact of anaerobic incubation and incubation with Broc on the composition of the gut microbiome.

Alpha Diversity

Fecal Culture Composition

We will start by plotting the composition of our samples before and after incubation with broccoli sprouts

Without Fecal Stocks
With Fecal Stocks
Animated

Alpha Diversity Metrics

Next let’s calculate alpha-diversity metrics and see if there are any significant differences following incubation. We will conduct a Kruskal-Wallis test followed by a Dunn’s Post-Hoc test.

adiv_dunn <- ps_raw %>% 
  subset_samples(treatment %in% c('fecal_stock', 'broc', 'NC')) %>%
  estimate_richness(measures = c('Observed','Shannon','Simpson')) %>%
  cbind(., data.frame(sample_data(ps_raw %>% 
                                    subset_samples(treatment %in% c('fecal_stock', 'broc', 'NC'))))) %>%
  pivot_longer(cols = c('Observed', 'Shannon', 'Simpson'), names_to = 'measure') %>%
  dplyr::select(measure, value, treatment, fecal_sample) %>%
  group_by(measure) %>%
  nest() 
walk2(.x = adiv_dunn$data, .y = list('Observed', 'Shannon', 'Simpson'), function(x, y) {
  print(y)
  dunn.test(x$value, g = x$treatment, kw = T, method = 'bh')})
## [1] "Observed"
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 19.4372, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## Col Mean-|
## Row Mean |       broc   fecal_st
## ---------+----------------------
## fecal_st |  -3.950143
##          |    0.0001*
##          |
##       NC |  -0.279431   3.670712
##          |     0.3900    0.0002*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
## [1] "Shannon"
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 19.3574, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## Col Mean-|
## Row Mean |       broc   fecal_st
## ---------+----------------------
## fecal_st |  -3.835403
##          |    0.0002*
##          |
##       NC |  -0.050800   3.784603
##          |     0.4797    0.0001*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
## [1] "Simpson"
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 19.3961, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## Col Mean-|
## Row Mean |       broc   fecal_st
## ---------+----------------------
## fecal_st |  -3.708403
##          |    0.0002*
##          |
##       NC |   0.203200   3.911603
##          |     0.4195    0.0001*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

For all measures of alpha diversity, there was a significant difference between the pre-incubation fecal stocks and both the Broc and NC samples. No significant differences in alpha diversity were found between the Broc and NC treatments, however, sample T4669 did change drastically in composition.

Beta Diversity

Next let’s analyze the beta diversity of our samples following incubation. We will start by calculating beta-diversity including the pre-incubation fecal stocks.

Beta Diversity with Fecal Stocks

ordBCall <- ordinate(ps_fbnv, method = 'PCoA', distance = 'bray')
allbd <- plot_ordination(ps_fbnv, ordBCall, color = 'treatment') +
  geom_point(aes(text = paste0('Fecal Sample: ', fecal_sample, '\n',
                               'Treatment: ', treatment)),
             size = 3) +
  theme_cowplot() +
  colTrt +
  ggtitle('Beta Diversity of All Samples')

ggplotly(allbd, tooltip = 'text', width = 900, height = 600)

We see three clusters forming: one cluster contains all the pre-incubation samples,and two additional clusters, each of which contained a mix of both Broc and NC samples.

k-Means Clustering

To verify the observed clusters, k-means clustering on the PCoA with 25 random starts and using between 2 and 10 clusters was conducted.

ugh <- ordinate(ps_fbnv, method = 'PCoA', distance = 'bray') %>%
  plot_ordination(ps_fbnv, ., justDF = TRUE) 


map_df(2:10, function(x) kmeans(ugh[,1:2], x, nstart = 25) %>% extract('tot.withinss')) %>%
  cbind(2:10, .) %>%
    plot(type = 'b', xlab = 'Number of Clusters K', ylab = 'Total Within-Clusters Sum of Squares')

The elbow method on within-clusters sum of squares verified the presence of 3 clusters within our bacterial samples. Next let’s replot our samples with the clustering:

ughc <- kmeans(ugh[,1:2], 3, nstart = 25) %>% 
  extract('cluster') %>%
  cbind(ugh, .) %>%
  modify_at('cluster', as.factor)


c <- ggplot(ughc, aes(x = Axis.1, y = Axis.2)) + 
  stat_ellipse(aes(group = cluster, fill = cluster), geom = 'polygon', color = 'black', alpha = 0.37, show.legend = FALSE) +
  #geom_point(aes(fill = treatment), color = 'black', pch = 21, size = 4, show.legend = TRUE) +
  xlab('Axis 1 [50.1%]') +
  ylab('Axis 2 [26.5%]') +
  theme_cowplot() + 
  theme(legend.position = 'none') +
  scale_x_continuous(limits = c(-0.8, 0.7)) +
  scale_y_continuous(limits = c(-0.8, 0.5)) +
  scale_fill_manual(values = c('#A50104', '#7DAF9C','#FCBA04')) +
  annotate('text', label = c('Fecal Stock Cluster', 'Mixed Cluster 1', 'Mixed Cluster 2'), 
           x = c(-0.5, 0.45, 0.1), y = c(0.3, 0.35, -0.8),
           size = 5)
  
trt <- RColorBrewer::brewer.pal(n = 3, name = 'Accent')

d <- ggplot(ughc, aes(x = Axis.1, y = Axis.2)) + 
  #stat_ellipse(aes(group = cluster, fill = cluster), geom = 'polygon', color = 'black', alpha = 0.37, show.legend = FALSE) +
  geom_point(aes(fill = treatment), color = 'black', pch = 21, size = 4, show.legend = TRUE) +
  xlab('Axis 1 [50.1%]') +
  ylab('Axis 2 [26.5%]') +
  theme_cowplot() + 
  scale_x_continuous(limits = c(-0.8, 0.7)) +
  scale_y_continuous(limits = c(-0.8, 0.5)) +
  scale_fill_manual(values = c(trt[2], trt[3], trt[1]), 
                    labels = c('Fecal Stock', 'NC', 'Broc'),
                    name = 'Treatment') +
  theme(legend.position = c(0.02,0.15))


aligned_plots <- align_plots(d, c, align = 'hv', axis = 'tbl')

ap <- ggdraw(aligned_plots[[2]]) + draw_plot(aligned_plots[[1]])

ap

As expected, one of these clusters was composed entirely of pre-incubation fecal stock samples, except for a single outlying NC sample. Consistent with PCoA findings, the other two clusters were composed of a mix of both NC and Broc samples. Next let’s remove the fecal stocks and replot with just Broc and NC.

Beta Diversity Without Fecal Stocks

s1 <- ordinate(ps_bnv, method = 'PCoA', distance = 'bray') %>%
  plot_ordination(ps_bnv, ., justDF = TRUE) 

map_df(1:9, function(x) kmeans(s1[,1:2], x, nstart = 25) %>% extract('tot.withinss')) %>%
  cbind(1:9, .) %>%
    plot(type = 'b', xlab = 'Number of Clusters K', ylab = 'Total Within-Clusters Sum of Squares')

s1c <- kmeans(s1[,1:2], 2, nstart = 25) %>% 
  extract('cluster') %>%
  cbind(s1, .) 
s1c$cluster %<>% as.factor()


sc <- ggplot(s1c, aes(x = Axis.1, y = Axis.2)) + 
  stat_ellipse(aes(group = cluster, fill = cluster), geom = 'polygon', color = 'black', alpha = 0.37, show.legend = FALSE) +
  #geom_point(aes(fill = treatment), color = 'black', pch = 21, size = 4, show.legend = TRUE) +
  xlab('Axis 1 [69.2%]') +
  ylab('Axis 2 [22.3%]') +
  theme_cowplot() + 
  theme(legend.position = 'none') +
  scale_x_continuous(limits = c(-0.5, 1)) +
  scale_y_continuous(limits = c(-1, 0.8)) +
  scale_fill_manual(values = rev(c('#A50104', '#7DAF9C'))) +
  annotate('text', label = c('Mixed Cluster 1', 'Mixed Cluster 2'), 
           x = c(-0.3, 0.3), y = c(0.25, -0.6),
           size = 5)

sd <- ggplot(s1c, aes(x = Axis.1, y = Axis.2)) + 
  #stat_ellipse(aes(group = cluster, fill = cluster), geom = 'polygon', color = 'black', alpha = 0.37, show.legend = FALSE) +
  geom_point(aes(fill = treatment), color = 'black', pch = 21, size = 4, show.legend = TRUE) +
  xlab('Axis 1 [69.2%]') +
  ylab('Axis 2 [22.3%]') +
  theme_cowplot() + 
  scale_x_continuous(limits = c(-0.5, 1)) +
  scale_y_continuous(limits = c(-1, 0.8)) +
  scale_fill_manual(values = c(trt[3], trt[1]), 
                    labels = c('NC', 'Broc'),
                    name = 'Treatment') +
  theme(legend.position = c(0.02,0.1))
sap <- align_plots(sd, sc, align = 'hv', axis = 'tbl')

ggdraw(sap[[2]]) + draw_plot(sap[[1]])

In the post-incubation PCoA, the two clusters previously observed persisted and were once again verified by k-means clustering. Next let’s see if the treatment we observe is due to treatment or fecal stock.

Treatment PERMANOVA

#First calculate the Bray-Curtis distances of the samples
distmat <- phyloseq::distance(ps_bnv, method = 'bray')
#Extract the metadata
mdata <- data.frame(sample_data(ps_bnv)) %>% dplyr::select(fecal_sample, treatment)
#Set the seed for reproducibility
set.seed(120)
#Run betadisper to verify the distrubtion of our groups is equal, an underlying assumption of PERMANOVA
bdisp <- betadisper(distmat, mdata$treatment)
#Evaluate using a permutation test
permutest(bdisp)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.01617 0.016169 0.3797    999   0.55
## Residuals 18 0.76655 0.042586
adonis2(distmat~treatment, data = mdata)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distmat ~ treatment, data = mdata)
##           Df SumOfSqs      R2      F Pr(>F)
## treatment  1   0.1741 0.04945 0.9365  0.366
## Residual  18   3.3456 0.95055              
## Total     19   3.5196 1.00000
plot(bdisp, label = F, col = paltrt[2:5]) %$%
legend('bottomright', legend = c('NC', 'broc'), col = paltrt[2:3], pch = 19) %$%
points(bdisp$centroids, pch = 20, cex = 2, col = paltrt[2:3])

Fecal Stock PERMANOVA

#First calculate the Bray-Curtis distances of the samples
distmat <- phyloseq::distance(ps_bnv, method = 'bray')
#Extract the metadata
mdata <- data.frame(sample_data(ps_bnv)) %>% dplyr::select(metabolomics_pos_sample, fecal_sample, treatment)
#Set the seed for reproducibility
set.seed(120)
#Run betadisper to verify the distrubtion of our groups is equal, an underlying assumption of PERMANOVA
bdisp <- betadisper(distmat, mdata$fecal_sample) #Significant but due PERMANOVA is robust to this when there is a balanced design (like what we have)

#Evaluate using a permutation test
permutest(bdisp)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq          F N.Perm Pr(>F)    
## Groups     9 0.19525 0.021694 1.9231e+30    999  0.001 ***
## Residuals 10 0.00000 0.000000                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(distmat~fecal_sample, data = mdata)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distmat ~ fecal_sample, data = mdata)
##              Df SumOfSqs      R2      F Pr(>F)    
## fecal_sample  9   2.8892 0.82088 5.0922  0.001 ***
## Residual     10   0.6304 0.17912                  
## Total        19   3.5196 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(bdisp, label = F, col = palsmp) %$%
legend('bottomright', legend = names(palsmp), col = palsmp, pch = 19) %$%
points(bdisp$centroids, pch = 20, cex = 2, col = palsmp)

Treatment was found to be non-significant by PERMANOVA, verifying our previous observation that the clustering observed was not driven by the incubation of broccoli with the fecal stocks. Conversely, sample donor was found to be significant by PERMANOVA. Cumulatively, these findings indicate that while anaerobic incubation significantly alters the composition of the microbiome, broccoli sprouts did not significantly alter the mi-crobiome relative to the negative control digest. Additionally, we identified two sub-populations within our samples, regardless of exposure to broccoli sprouts.

Gut Microbiome Composition Influences Nitrile Production

Using mass spectrometry, we measured the abundance of SFN-NIT and IBN-NIT in our samples.

Nitriles Present in Each Cluster

mdataveg <- data.frame(sample_data(ps_broc)) %>% dplyr::select(metabolomics_pos_sample, fecal_sample, treatment)
metabinterest <- read_csv('./data/metab_interest.csv') %>%
  right_join(., mdataveg, by = c('sample' = 'metabolomics_pos_sample'))

nitdata <- metabinterest %>% #Quantified IBN, IBNNIT, and SFNNIT from multiquant
  pivot_longer(where(is.numeric), names_to = 'metabolite', values_to = 'intensity') %>%
  left_join(., data.frame(sample_data(ps_broc)) %>% select(fecal_sample, group)) %>%
  filter(treatment == 'broc')

kt <- nitdata %>%
  group_by(metabolite) %>%
  nest() %>%
  mutate('p-value' = map_dbl(data, function(x) kruskal.test(intensity ~ group, data = x) %>% broom::tidy() %>% pull(p.value))) %>%
  dplyr::select(-data) %>%
  rename('Metabolite' = metabolite)
  
knitr::kable(kt, caption = 'Differences in Means NIT between Clusters')
Differences in Means NIT between Clusters
Metabolite p-value
IBNNIT 0.0090234
SFNNIT 0.0090234 
SEM <- function(x) sd(x)/sqrt(length(x))

nit2 <- nitdata %>%
  group_by(metabolite, group) %>%
  summarise(mean = mean(intensity), se = SEM(intensity)) %>%
  modify_at('metabolite', as.factor) %>%
  modify_at('metabolite', fct_relevel, c('SFNNIT', 'IBNNIT'))

sfnnit <- nit2 %>%
  filter(metabolite == 'SFNNIT') %>%
  modify_at('group', function(x) as.factor(ifelse(x == 'A', 1, 2)))

ibnnit <- nit2 %>%
  filter(metabolite == 'IBNNIT') %>%
  modify_at('group', function(x) as.factor(ifelse(x == 'A', 1, 2)))

s <- ggplot(sfnnit, aes(x = group, y = mean, group = group)) +
  geom_errorbar(aes(ymax = mean + se, ymin = mean - se), width = 0.5) +
  geom_col(fill = 'darkgrey', color = 'black') +
  ylab('SFN-NIT (Relative Intensity)') +
  xlab('Cluster') +
  theme_cowplot() +
  scale_fill_manual(labels = c('Cluster 1', 'Cluster 2'),
                    name = 'Cluster') +
  scale_y_continuous(label = scales::label_comma(),
                     expand = c(0,0),
                     limits = c(0,172000)) +
  annotate('text', label = c('*'), 
           x = 2, y = 35000,
           size = 10) 

i <- ggplot(ibnnit, aes(x = group, y = mean, group = group)) +
  geom_errorbar(aes(ymax = mean + se, ymin = mean - se), width = 0.5) +
  geom_col(fill = 'darkgrey', color = 'black') +
  ylab('IBN-NIT (Relative Intensity)') +
  xlab('Cluster') +
  theme_cowplot() +
  scale_fill_manual(labels = c('Cluster 1', 'Cluster 2'),
                    name = 'Cluster') +
  scale_y_continuous(label = scales::label_comma(),
                     expand = c(0,0),
                     limits = c(0,13000)) +
  annotate('text', label = c('*'), 
           x = 2, y = 6300,
           size = 10) 

plot_grid(s, i, labels = c('A','B'), axis = 'b')

We find that the relative abundance of SFN-NIT in fecal culture medium, measured by integrating their chromatographic peak areas, was significantly different (p = 0.009) between samples belonging to clusters 1 and 2, with the relative intensity being 156,400 and 15,336 for each cluster, respectively. It was also significantly different (p = 0.009) for IBN-NIT, with the relative intensity being 11,844 and 4,399 for each cluster, respectively. Next we will project the abundance of each NIT onto our beta-diversity plot and run a PERMANOVA.

Nitriles PERMANOVA

ordBCveg <- ordinate(ps_broc, method = 'PCoA', distance = 'bray') %>%
  plot_ordination(ps_broc, ., justDF = TRUE) %>%
  left_join(., metabinterest, by = c('metabolomics_pos_sample' = 'sample')) %>%
  select(-ends_with('.y'))


vegbd <- ggplot(ordBCveg, aes(x = Axis.1, y = Axis.2, size = SFNNIT, color = IBNNIT)) +
  geom_point(aes(text = paste0('Fecal Sample: ', fecal_sample.x, '\n',
                               'Treatment: ', treatment.x, '\n',
                               'IBNNIT: ', IBNNIT, '\n',
                               'SFNNIT: ', SFNNIT))) +
  scale_color_viridis_c() +
  theme_cowplot() +
  ggtitle('Beta Diversity of Digests Following Incubation')

Interactive Plot
Static Plot

set.seed(101)
mdataveg <- data.frame(sample_data(ps_broc)) %>% dplyr::select(metabolomics_pos_sample, fecal_sample, treatment)
distmatveg <- phyloseq::distance(ps_broc, method = 'bray')
metabinterest <- read_csv('./data/metab_interest.csv') %>%
  right_join(., mdataveg, by = c('sample' = 'metabolomics_pos_sample'))

suppressWarnings(    #SuppressWarnings because broom::tidy() does not play well PERMANOVA outputs
perma <- metabinterest %>%
  pivot_longer(where(is.numeric), names_to = 'Metabolite', values_to = 'intensity') %>%
  dplyr::group_by(Metabolite) %>%
  nest() %>%
  mutate(permanova = map(data, function(x){ 
    adonis2(distmatveg~intensity, data = x) %>% 
      broom::tidy() %>%
      magrittr::extract('p.value') %>%
      drop_na() 
    })) %>%
  unnest(cols = 'permanova') %>% 
  select(-data)
)

knitr::kable(perma, caption = 'PERMANOVA Results for NITs')
PERMANOVA Results for NITs
Metabolite p.value
IBNNIT 0.018
SFNNIT 0.008

PERMANOVA on the broccoli samples alone further corroborated that both SFN-NIT and IBN-NIT were significantly different between the two clusters, indicating that differing microbiome compositions yielded differing metabolic capabilities. These results indicate that the composition of the gut microbiome influences con-version of GLS to NITs.

Members of Family Clostridiaceae Are Associated with Nitrile

LEfSe Analysis

To identify which taxa discriminate between the two observed clusters (High vs Low-NIT producers), we conducted a Linear Discriminant Analysis Effect Size (LEfSe). LEfSe was run on the browser-based Galaxy Module.

LDA Results

Cladogram

In the post-incubation samples, LEfSe identified the family Enterobacteriaceae and the order Bacteroidales as being associated with the low-NIT producing (cluster 2) samples. Conversely, the phylum Firmicutes, specifically members of the Clostridia order and Clostridiaceae family, were associated with high-NIT producing (cluster 1) individuals.

Negative Binomial Model

To verify these results and identify individual taxa differentially abundant between high- and low-NIT producers, a negative binomial model was built for each taxa with abundance as the response variable and the interaction of treatment and group as the predictors.

library(lme4)
library(phia)
set.seed(123)

psglm <- ps_counts %>%
  subset_samples(treatment %in% c('broc', 'NC')) %>%
  rarefy_even_depth() %>%
  psmelt()

glmnb <- psglm %>%
  group_by(OTU) %>%
  nest() %>%
  mutate_at( # change the data in 'lmer' column
    "data",
    purrr::map, # allows iteration over data in 'lmer' column, returns list
    function(x) {
      lmer <- try(glmer.nb(Abundance ~ treatment*group + (1 |fecal_sample), data = x), silent = TRUE) # do the lmer test
      c(lmer, x) # return results of lmer, as well as data for (r)anova later
    } 
  )

#Remove the failed runs which contain error messages
glmnbclean <- glmnb[sapply(glmnb$data, function(x) class(x[[1]]) !=  'character'),]

#Function to calculate the log2FoldChange of each microbial species:
#Helper function
l2fcmin <- function(x, counts, factor, class1, class2){
  x[[counts]] <- log2(x[[counts]]+sqrt(nzmin(x[[counts]])*0.001))
  c1 <- x[[counts]][sapply(x[[factor]], function(x) x == class1)]
  c2 <- x[[counts]][sapply(x[[factor]], function(x) x == class2)]
  lf2 <- mean(c1[sapply(c1,is.finite)]) - mean(c2[sapply(c2, is.finite)])
  return(lf2)
}

nzmin <- function(x){
  min(x[x > 0])
}


lfc <- psglm %>%
  group_by(OTU, treatment) %>%
  nest() %>%
  mutate(L2FC = purrr::map(data, function(x) l2fcmin(x, counts = 'Abundance', factor = 'group', class1 = 'A', class2 = 'B'))) %>%
  dplyr::select(-data) %>%
  unnest(L2FC) %>%
  ungroup()

glmres <- glmnbclean %>%
  group_by(OTU) %>%
  group_modify(~{
    testInteractions(.x$data[[1]][[1]], fixed = 'treatment', across = 'group') %>% broom::tidy()
  }) %>%
  ungroup() %>%
  mutate(padj = p.adjust(p.value, method = 'BH')) %>%
  left_join(lfc, by = c('OTU' = 'OTU', 'term' = 'treatment'))

glmsig <- glmres %>%
  filter(padj <= 0.05) %>%
  left_join(data.frame(tax_table(ps_counts)) %>% rownames_to_column('OTU')) 
datatable(glmsig,
           filter = 'top',
           extensions = 'Buttons',
           options = list(pageLength = 8,
                          dom = 'Bfrtip',
                          scrollX = TRUE))

7 taxa were significantly different (p < 0.05) between the observed clusters all of which belonged to the families Enterobacteriaceae, Enterococcaceae, Clostridiaceae, and Cori-obacteriaceae (Table 1). Members of Enterobacteriaceae were significantly more abundant in low-NIT individuals while taxa from Clostridiaceae and Coriobacteriaceae were more abundant in high-NIT individuals. Members of family Enterococcaceae were found to be more abundant in low NIT producers in NC samples but more abundant in high NIT producers.

Correlation Analysis

Lastly, spearman’s correlations between the abundance of each family and the abundance of NIT across all Broc samples were performed.

psfam <- ps_broc %>%
  tax_glom(taxrank = 'Family') %>%
  psmelt()

cordata <- psfam %>%
  dplyr::rename('ASV' = 'OTU') %>%
  dplyr::select(ASV, Abundance, metabolomics_pos_sample, group, fecal_sample, Family) %>%
  dplyr::rename('sample' = 'metabolomics_pos_sample') %>%
  left_join(metabinterest) %>%
  pivot_longer(cols = c('IBNNIT', 'SFNNIT'), names_to = 'metabolite', values_to = 'intensity') %>%
  group_by(Family, metabolite) %>%
  nest() %>%
  mutate(rho = map(data, function(x){ 
    stats::cor(x$intensity, x$Abundance, method = 'spearman')
    })) %>%
  unnest(cols = 'rho') %>%
  ungroup() %>%
  dplyr::select(-data) %>%
  drop_na()

hmap <- ggplot(cordata, aes(x = Family, y = metabolite, fill = rho)) +
  geom_tile(aes(text = paste0('Metabolite: ', metabolite, '\n',
                                'Family: ', Family, '\n',
                                'Rho:',  round(rho,3)))) +
  coord_fixed() +
  scale_fill_gradientn(colors = c('blue', 'white', 'red')) +
  theme_cowplot() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_x_discrete(expand = c(0,0)) +
  scale_y_discrete(expand = c(0,0)) +
  ggtitle('Spearman Correlations Between Families and Metabolites for Broccoli')

ggplotly(hmap, tooltip = 'text', height = 450, width = 900)

Enterobacteriaceae was found to be strongly negatively correlated with both SFN-NIT and IBN-NIT while the Clostridiaceae family was strongly positively correlated with SFN-NIT and IBN-NIT. Enterococcaceae had a moderately strong correlation with both SFN-NIT and IBN-NIT.

SFN-NIT Variation Exists in Humans

Data from our fecal culture experiments and results from others suggest that the microbiome may have an important role in how glucoraphanin is metabolized to SFN-NIT but to date it is not clear how much SFN-NIT is excreted following broccoli sprout con-sumption by humans. To determine if the variations in SFN-NIT production we observed in fecal cultures are also observed with human broccoli consumption, we analyzed urine samples from participants fed a single dose of broccoli sprouts containing 200 µmol of SFN equivalents (full results published in [16]).

ad <- rawdata_atwell%>%
  dplyr::select(-phase, -sample_site) %>%
  group_by(subject, time) %>%
  mutate(SFN_tot = SFN+SFNNIT+SFN_CYS+SFN_NAC+SFN_CG+SFN_GSH) %>%
  modify_at(c('time','subject'), as.factor)

v <- volumes %>%
  #Convert mL to L
  mutate(across(c(2:7), ~.x/1000)) %>%
  pivot_longer(cols = -subject, names_to = 'time', values_to = 'volume') %>%
  modify_at('time', as.integer) %>%
  modify_at(c('time','subject'), as.factor) %>%
  #Add µM to volume data
  left_join(ad) %>%
  group_by(subject, time) %>%
  #Multiple each µM by volume of urine to get µmol
  mutate(across(starts_with('SFN'), ~.x*volume)) 

bleh <- v %>%
  dplyr::select(subject, time, SFNNIT, SFN_NAC) %>%
  pivot_longer(where(is.numeric), names_to = 'metabolite', values_to = 'umol')

test <- ggplot(bleh, aes(x = time, y = umol, color = subject, group = subject)) +
  geom_bar(aes(group = time),stat = 'summary', fun = mean, fill = 'grey', color = 'black', alpha = 0.5) +
  geom_point(size = 3, aes(text = paste0('Metabolite: ', metabolite, '\n',
                                         'µmol: ', umol, '\n',
                                         'Subject: ', subject, '\n',
                                         'Time: ', time, 'h'))) +
  geom_path(size = 1.5) + 
  theme_minimal_hgrid() +
  ylab('SFN-NAC (µmol)') +
  xlab('Time (Hours)') + 
  scale_y_continuous(limits = c(0, 65), expand = c(0.00,0.00), breaks = seq(0, 70, by = 5)) +
  scale_x_discrete(expand = c(0.005,0.005)) +
  scale_color_manual(name = 'Subject ID', 
                     values = c('#F8766D', '#D89000', '#A3A500', '#39B600', '#00BF7D', '#00BFC4', '#00B0F6', '#9590FF', '#E76BF3', '#FF62BC')) +
  facet_wrap(~metabolite)

ggplotly(test, tooltip = 'text')

Conclusions:

Human clinical trials utilizing broccoli sprouts, and other whole food interventions, frequently observe high levels of inter-individual variation between subjects coinciding with equivocal efficacy. The source of variation is still unknown, but a growing body of evidence points towards the gut microbiome as playing a key role in dictating the metabolic fate of GLS from whole foods. The role of the gut microbiome on GLS metab-olism using dynamic and complex systems is still relatively undescribed and the human excretion of NITs, a metabolite considered to be biologically inert of GLS, from broccoli sprouts has yet to be measured. Utilizing urine samples from humans fed 200 µmol of SFN-equivalents in broccoli sprouts, we detected high levels of inter-individual variation in SFN-NIT excretion following consumption. Using an ex vivo fecal incubation system, we determined that the gut microbiome influences production of NITs from broccoli sprouts. Specifically, we found that individuals with gut microbiomes enriched with members of family Clostridiaceae had relatively high levels of SFN-NIT and IBN-NIT while those with gut microbiomes abundant with taxa from the family Enterobacteriaceae had a relatively low abundance of NITs. Cumulatively, these findings point towards the gut microbiome as playing a key role in influencing the production of NITs from broccoli sprouts. Un-derstanding factors which drive variation in GLS metabolism may explain why some individuals benefit more greatly from the anti-cancer properties of cruciferous vegetables than others. For a full discussion of our results please see our article in Nutrients.

References

  1. Holst, B.; Williamson, G. A Critical Review of the Bioavailability of Glucosinolates and Related Compounds. Nat. Prod. Rep. 2004, 21, 425–447, doi:10.1039/B204039P.
  2. Verkerk, R.; Schreiner, M.; Krumbein, A.; Ciska, E.; Holst, B.; Rowland, I.; Schrijver, R.D.; Hansen, M.; Gerhäuser, C.; Mithen, R.; et al. Glucosinolates in Brassica Vegetables: The Influence of the Food Supply Chain on Intake, Bioavailability and Human Health. Molecular Nutrition & Food Research 2009, 53, S219–S219, doi:10.1002/mnfr.200800065.
  3. Abbaoui, B.; Lucas, C.R.; Riedl, K.M.; Clinton, S.K.; Mortazavi, A. Cruciferous Vegetables, Isothiocyanates and Bladder Cancer Prevention. Mol Nutr Food Res 2018, 62, e1800079, doi:10.1002/mnfr.201800079.
  4. Ambrosone, C.B.; McCann, S.E.; Freudenheim, J.L.; Marshall, J.R.; Zhang, Y.; Shields, P.G. Breast Cancer Risk in Premenopausal Women Is Inversely Associated with Consumption of Broccoli, a Source of Isothiocyanates, but Is Not Modified by GST Genotype. J Nutr 2004, 134, 1134–1138, doi:10.1093/jn/134.5.1134.
  5. Aune, D.; Giovannucci, E.; Boffetta, P.; Fadnes, L.T.; Keum, N.; Norat, T.; Greenwood, D.C.; Riboli, E.; Vatten, L.J.; Tonstad, S. Fruit and Vegetable Intake and the Risk of Cardiovascular Disease, Total Cancer and All-Cause Mortality—a Systematic Review and Dose-Response Meta-Analysis of Prospective Studies. Int J Epidemiol 2017, 46, 1029–1056, doi:10.1093/ije/dyw319.
  6. Lam, T.K.; Gallicchio, L.; Boyd, K.; Shiels, M.; Hammond, E.; Tao, X. (Grant); Chen, L.; Robinson, K.A.; Caulfield, L.E.; Herman, J.G.; et al. Cruciferous Vegetable Consumption and Lung Cancer Risk: A Systematic Review. Cancer Epidemiol Biomarkers Prev 2009, 18, 184–195, doi:10.1158/1055-9965.EPI-08-0710.
  7. Basten, G.P.; Bao, Y.; Williamson, G. Sulforaphane and Its Glutathione Conjugate but Not Sulforaphane Nitrile Induce UDP-Glucuronosyl Transferase (UGT1A1) and Glutathione Transferase (GSTA1) in Cultured Cells. Carcinogenesis 2002, 23, 1399–1404, doi:10.1093/carcin/23.8.1399.
  8. Azarenko, O.; Jordan, M.A.; Wilson, L. Erucin, the Major Isothiocyanate in Arugula (Eruca Sativa), Inhibits Proliferation of MCF7 Tumor Cells by Suppressing Microtubule Dynamics. PLoS One 2014, 9, doi:10.1371/journal.pone.0100599.
  9. Hanlon, N.; Coldham, N.; Sauer, M.J.; Ioannides, C. Up-Regulation of the CYP1 Family in Rat and Human Liver by the Aliphatic Isothiocyanates Erucin and Sulforaphane. Toxicology 2008, 252, 92–98, doi:10.1016/j.tox.2008.08.002.
  10. Jakubíková, J.; Sedlák, J.; Mithen, R.; Bao, Y. Role of PI3K/Akt and MEK/ERK Signaling Pathways in Sulforaphane- and Erucin-Induced Phase II Enzymes and MRP2 Transcription, G2/M Arrest and Cell Death in Caco-2 Cells. Biochemical Pharmacology 2005, 69, 1543–1552, doi:10.1016/j.bcp.2005.03.015.
  11. Zhang, Y.; Talalay, P.; Cho, C.G.; Posner, G.H. A Major Inducer of Anticarcinogenic Protective Enzymes from Broccoli: Isolation and Elucidation of Structure. Proc Natl Acad Sci U S A 1992, 89, 2399–2403.
  12. McWalter, G.K.; Higgins, L.G.; McLellan, L.I.; Henderson, C.J.; Song, L.; Thornalley, P.J.; Itoh, K.; Yamamoto, M.; Hayes, J.D. Transcription Factor Nrf2 Is Essential for Induction of NAD(P)H:Quinone Oxidoreductase 1, Glutathione S-Transferases, and Glutamate Cysteine Ligase by Broccoli Seeds and Isothiocyanates. J Nutr 2004, 134, 3499S-3506S, doi:10.1093/jn/134.12.3499S.
  13. Bonnesen, C.; Eggleston, I.M.; Hayes, J.D. Dietary Indoles and Isothiocyanates That Are Generated from Cruciferous Vegetables Can Both Stimulate Apoptosis and Confer Protection against DNA Damage in Human Colon Cell Lines. Cancer Res 2001, 61, 6120–6130.
  14. Clarke, J.D.; Riedl, K.; Bella, D.; Schwartz, S.J.; Stevens, J.F.; Ho, E. Comparison of Isothiocyanate Metabolite Levels and Histone Deacetylase Activity in Human Subjects Consuming Broccoli Sprouts or Broccoli Supplement. J Agric Food Chem 2011, 59, 10955–10963, doi:10.1021/jf202887c.
  15. Cao, C.; Wu, H.; Vasilatos, S.N.; Chandran, U.; Qin, Y.; Wan, Y.; Oesterreich, S.; Davidson, N.E.; Huang, Y. HDAC5-LSD1 Axis Regulates Antineoplastic Effect of Natural HDAC Inhibitor Sulforaphane in Human Breast Cancer Cells. Int J Cancer 2018, 143, 1388–1401, doi:10.1002/ijc.31419.
  16. Atwell, L.L.; Hsu, A.; Wong, C.P.; Stevens, J.F.; Bella, D.; Yu, T.-W.; Pereira, C.B.; Löhr, C.V.; Christensen, J.M.; Dashwood, R.H.; et al. Absorption and Chemopreventive Targets of Sulforaphane in Humans Following Consumption of Broccoli Sprouts or a Myrosinase-Treated Broccoli Sprout Extract. Mol Nutr Food Res 2015, 59, 424–433, doi:10.1002/mnfr.201400674.
  17. Yagishita, Y.; Fahey, J.W.; Dinkova-Kostova, A.T.; Kensler, T.W. Broccoli or Sulforaphane: Is It the Source or Dose That Matters? Molecules 2019, 24, doi:10.3390/molecules24193593.
  18. Epplein, M.; Wilkens, L.R.; Tiirikainen, M.; Dyba, M.; Chung, F.-L.; Goodman, M.T.; Murphy, S.P.; Henderson, B.E.; Kolonel, L.N.; Marchand, L.L. Urinary Isothiocyanates; Glutathione S-Transferase M1, T1, and P1 Polymorphisms; and Risk of Colorectal Cancer: The Multiethnic Cohort Study. Cancer Epidemiol Biomarkers Prev 2009, 18, 314–320, doi:10.1158/1055-9965.EPI-08-0627.
  19. Lin, H.J.; Probst-Hensch, N.M.; Louie, A.D.; Kau, I.H.; Witte, J.S.; Ingles, S.A.; Frankl, H.D.; Lee, E.R.; Haile, R.W. Glutathione Transferase Null Genotype, Broccoli, and Lower Prevalence of Colorectal Adenomas. Cancer Epidemiol Biomarkers Prev 1998, 7, 647–652.
  20. Moy, K.A.; Yuan, J.-M.; Chung, F.-L.; Van Den Berg, D.; Wang, R.; Gao, Y.-T.; Yu, M.C. Urinary Total Isothiocyanates and Colorectal Cancer: A Prospective Study of Men in Shanghai, China. Cancer Epidemiol Biomarkers Prev 2008, 17, 1354–1359, doi:10.1158/1055-9965.EPI-07-2841.
  21. Slattery, M.L.; Kampman, E.; Samowitz, W.; Caan, B.J.; Potter, J.D. Interplay between Dietary Inducers of GST and the GSTM-1 Genotype in Colon Cancer. International Journal of Cancer 2000, 87, 728–733, doi:10.1002/1097-0215(20000901)87:5<728::AID-IJC16>3.0.CO;2-G.
  22. Seow, A.; Vainio, H.; Yu, M.C. Effect of Glutathione-S-Transferase Polymorphisms on the Cancer Preventive Potential of Isothiocyanates: An Epidemiological Perspective. Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis 2005, 592, 58–67, doi:10.1016/j.mrfmmm.2005.06.004.
  23. Seow, A.; Yuan, J.-M.; Sun, C.-L.; Van Den Berg, D.; Lee, H.-P.; Yu, M.C. Dietary Isothiocyanates, Glutathione S -Transferase Polymorphisms and Colorectal Cancer Risk in the Singapore Chinese Health Study. Carcinogenesis 2002, 23, 2055–2061, doi:10.1093/carcin/23.12.2055.
  24. Turner, F.; Smith, G.; Sachse, C.; Lightfoot, T.; Garner, R.C.; Wolf, C.R.; Forman, D.; Bishop, D.T.; Barrett, J.H. Vegetable, Fruit and Meat Consumption and Potential Risk Modifying Genes in Relation to Colorectal Cancer. International Journal of Cancer 2004, 112, 259–264, doi:10.1002/ijc.20404.
  25. Tijhuis, M.J.; Wark, P.A.; Aarts, J.M.M.J.G.; Visker, M.H.P.W.; Nagengast, F.M.; Kok, F.J.; Kampman, E. GSTP1 and GSTA1 Polymorphisms Interact with Cruciferous Vegetable Intake in Colorectal Adenoma Risk. Cancer Epidemiol Biomarkers Prev 2005, 14, 2943–2951.
  26. Yang, G.; Gao, Y.-T.; Shu, X.-O.; Cai, Q.; Li, G.-L.; Li, H.-L.; Ji, B.-T.; Rothman, N.; Dyba, M.; Xiang, Y.-B.; et al. Isothiocyanate Exposure, Glutathione S-Transferase Polymorphisms, and Colorectal Cancer Risk1234. Am J Clin Nutr 2010, 91, 704–711, doi:10.3945/ajcn.2009.28683.
  27. Brauer, H.A.; Libby, T.E.; Mitchell, B.L.; Li, L.; Chen, C.; Randolph, T.W.; Yasui, Y.Y.; Lampe, J.W.; Lampe, P.D. Cruciferous Vegetable Supplementation in a Controlled Diet Study Alters the Serum Peptidome in a GSTM1-Genotype Dependent Manner. Nutr J 2011, 10, 11, doi:10.1186/1475-2891-10-11.
  28. Navarro, S.L.; Chang, J.-L.; Peterson, S.; Chen, C.; King, I.B.; Schwarz, Y.; Li, S.S.; Li, L.; Potter, J.D.; Lampe, J.W. Modulation of Human Serum Glutathione S-Transferase-A1/2 Concentration by Cruciferous Vegetables in a Controlled Feeding Study Is Influenced by GSTM1 and GSTT1 Genotypes. Cancer Epidemiol Biomarkers Prev 2009, 18, 2974–2978, doi:10.1158/1055-9965.EPI-09-0701.
  29. Navarro, S.L.; Schwarz, Y.; Song, X.; Wang, C.-Y.; Chen, C.; Trudo, S.P.; Kristal, A.R.; Kratz, M.; Eaton, D.L.; Lampe, J.W. Cruciferous Vegetables Have Variable Effects on Biomarkers of Systemic Inflammation in a Randomized Controlled Trial in Healthy Young Adults12. J Nutr 2014, 144, 1850–1857, doi:10.3945/jn.114.197434.
  30. Peterson, S.; Schwarz, Y.; Li, S.S.; Li, L.; King, I.B.; Chen, C.; Eaton, D.L.; Potter, J.D.; Lampe, J.W. CYP1A2, GSTM1, and GSTT1 Polymorphisms and Diet Effects on CYP1A2 Activity in a Crossover Feeding Trial. Cancer Epidemiol Biomarkers Prev 2009, 18, 3118–3125, doi:10.1158/1055-9965.EPI-09-0589.
  31. Clarke, J.D.; Hsu, A.; Williams, D.E.; Dashwood, R.H.; Stevens, J.F.; Ho, E. Metabolism and Tissue Distribution of Sulforaphane in Nrf2 Knockout and Wild-Type Mice. Pharm Res 2011, 28, 3171–3179, doi:10.1007/s11095-011-0500-z.
  32. Clarke, J.D.; Riedl, K.; Bella, D.; Schwartz, S.J.; Stevens, J.F.; Ho, E. Comparison of Isothiocyanate Metabolite Levels and Histone Deacetylase Activity in Human Subjects Consuming Broccoli Sprouts or Broccoli Supplement. J Agric Food Chem 2011, 59, 10955–10963, doi:10.1021/jf202887c.
  33. Saha, S.; Hollands, W.; Teucher, B.; Needs, P.W.; Narbad, A.; Ortori, C.A.; Barrett, D.A.; Rossiter, J.T.; Mithen, R.F.; Kroon, P.A. Isothiocyanate Concentrations and Interconversion of Sulforaphane to Erucin in Human Subjects after Consumption of Commercial Frozen Broccoli Compared to Fresh Broccoli. Molecular Nutrition & Food Research 2012, 56, 1906–1916, doi:10.1002/mnfr.201200225.
  34. Keck, A.-S.; Finley, J.W. Cruciferous Vegetables: Cancer Protective Mechanisms of Glucosinolate Hydrolysis Products and Selenium. Integr Cancer Ther 2004, 3, 5–12, doi:10.1177/1534735403261831.
  35. Getahun, S.M.; Chung, F.-L. Conversion of Glucosinolates to Isothiocyanates in Humans after Ingestion of Cooked Watercress. Cancer Epidemiol Biomarkers Prev 1999, 8, 447–451.
  36. Li, F.; Hullar, M.A.J.; Schwarz, Y.; Lampe, J.W. Human Gut Bacterial Communities Are Altered by Addition of Cruciferous Vegetables to a Controlled Fruit- and Vegetable-Free Diet. J Nutr 2009, 139, 1685–1691, doi:10.3945/jn.109.108191.
  37. Li, F.; Hullar, M.A.J.; Beresford, S.A.A.; Lampe, J.W. Variation of Glucoraphanin Metabolism in Vivo and Ex Vivo by Human Gut Bacteria. British Journal of Nutrition 2011, 106, 408–416, doi:10.1017/S0007114511000274.
  38. Mullaney, J.A.; Kelly, W.J.; McGhie, T.K.; Ansell, J.; Heyes, J.A. Lactic Acid Bacteria Convert Glucosinolates to Nitriles Efficiently Yet Differently from Enterobacteriaceae. J. Agric. Food Chem. 2013, 61, 3039–3046, doi:10.1021/jf305442j.
  39. Luang‐In, V.; Narbad, A.; Nueno‐Palop, C.; Mithen, R.; Bennett, M.; Rossiter, J.T. The Metabolism of Methylsulfinylalkyl- and Methylthioalkyl-Glucosinolates by a Selection of Human Gut Bacteria. Molecular Nutrition & Food Research 2014, 58, 875–883, doi:10.1002/mnfr.201300377.
  40. Luang-In, V.; Albaser, A.A.; Nueno-Palop, C.; Bennett, M.H.; Narbad, A.; Rossiter, J.T. Glucosinolate and Desulfo-Glucosinolate Metabolism by a Selection of Human Gut Bacteria. Curr Microbiol 2016, 73, 442–451, doi:10.1007/s00284-016-1079-8.
  41. Luang-In, V.; Deeseenthum, S.; Udomwong, P.; Saengha, W.; Gregori, M. Formation of Sulforaphane and Iberin Products from Thai Cabbage Fermented by Myrosinase-Positive Bacteria. Molecules 2018, 23, 955, doi:10.3390/molecules23040955.
  42. Liu, X.; Wang, Y.; Hoeflinger, J.L.; Neme, B.P.; Jeffery, E.H.; Miller, M.J. Dietary Broccoli Alters Rat Cecal Microbiota to Improve Glucoraphanin Hydrolysis to Bioactive Isothiocyanates. Nutrients 2017, 9, doi:10.3390/nu9030262.
  43. Center for Food Safety and Applied Nutrition Guidance for Industry: Compliance with and Recommendations for Implementation of the Standards for the Growing, Harvesting, Packing, and Holding of Produce for Human Consumption for Sprout Operations 2017.
  44. Microbiological Safety Evaluations and Recommendations on Sprouted Seeds. International Journal of Food Microbiology 1999, 52, 123–153, doi:10.1016/S0168-1605(99)00135-X.
  45. Fahey, J.W.; Ourisson, P.J.; Degnan, F.H. Pathogen Detection, Testing, and Control in Fresh Broccoli Sprouts. Nutr J 2006, 5, 13, doi:10.1186/1475-2891-5-13.
  46. Henson, W.Y. U.S. EPA, Pesticides, Label, ECR CALCIUM HYPOCHLORITE T, 3/17/2011 2011.
  47. Aura, A.-M.; Härkönen, H.; Fabritius, M.; Poutanen, K. Development of AnIn VitroEnzymic Digestion Method for Removal of Starch and Protein and Assessment of Its Performance Using Rye AndWheat Breads. Journal of Cereal Science 1999, 29, 139–152, doi:10.1006/jcrs.1998.0229.
  48. Gil-Izquierdo, A.; Zafrilla, P.; Tomás-Barberán, F.A. An in Vitro Method to Simulate Phenolic Compound Release from the Food Matrix in the Gastrointestinal Tract. Eur Food Res Technol 2002, 214, 155–159, doi:10.1007/s00217-001-0428-3.
  49. Vallejo, F.; Gil-Izquierdo, A.; Pérez-Vicente, A.; García-Viguera, C. In Vitro Gastrointestinal Digestion Study of Broccoli Inflorescence Phenolic Compounds, Glucosinolates, and Vitamin C. J Agric Food Chem 2004, 52, 135–138, doi:10.1021/jf0305128.
  50. Sarvan, I.; Kramer, E.; Bouwmeester, H.; Dekker, M.; Verkerk, R. Sulforaphane Formation and Bioaccessibility Are More Affected by Steaming Time than Meal Composition during in Vitro Digestion of Broccoli. Food Chem 2017, 214, 580–586, doi:10.1016/j.foodchem.2016.07.111.
  51. Rychlik, J.; Olejnik, A.; Olkowicz, M.; Kowalska, K.; Juzwa, W.; Myszka, K.; Dembczyński, R.; Moyer, M.P.; Grajek, W. Antioxidant Capacity of Broccoli Sprouts Subjected to Gastrointestinal Digestion. J Sci Food Agric 2015, 95, 1892–1902, doi:10.1002/jsfa.6895.
  52. Guadamuro, L.; Dohrmann, A.B.; Tebbe, C.C.; Mayo, B.; Delgado, S. Bacterial Communities and Metabolic Activity of Faecal Cultures from Equol Producer and Non-Producer Menopausal Women under Treatment with Soy Isoflavones. BMC Microbiol 2017, 17, 93, doi:10.1186/s12866-017-1001-y.
  53. Caporaso, J.G.; Lauber, C.L.; Walters, W.A.; Berg-Lyons, D.; Huntley, J.; Fierer, N.; Owens, S.M.; Betley, J.; Fraser, L.; Bauer, M.; et al. Ultra-High-Throughput Microbial Community Analysis on the Illumina HiSeq and MiSeq Platforms. ISME J 2012, 6, 1621–1624, doi:10.1038/ismej.2012.8.
  54. Illumina 16S Sample Preparation Guide.
  55. Callahan, B.J.; McMurdie, P.J.; Rosen, M.J.; Han, A.W.; Johnson, A.J.A.; Holmes, S.P. DADA2: High-Resolution Sample Inference from Illumina Amplicon Data. Nature Methods 2016, 13, 581–583, doi:10.1038/nmeth.3869.
  56. Quast, C.; Pruesse, E.; Yilmaz, P.; Gerken, J.; Schweer, T.; Yarza, P.; Peplies, J.; Glöckner, F.O. The SILVA Ribosomal RNA Gene Database Project: Improved Data Processing and Web-Based Tools. Nucleic Acids Res 2013, 41, D590–D596, doi:10.1093/nar/gks1219.
  57. Axton, E.R.; Beaver, L.M.; St Mary, L.; Truong, L.; Logan, C.R.; Spagnoli, S.; Prater, M.C.; Keller, R.M.; Garcia-Jaramillo, M.; Ehrlicher, S.E.; et al. Treatment with Nitrate, but Not Nitrite, Lowers the Oxygen Cost of Exercise and Decreases Glycolytic Intermediates While Increasing Fatty Acid Metabolites in Exercised Zebrafish. J Nutr 2019, 149, 2120–2132, doi:10.1093/jn/nxz202.
  58. Housley, L.; Magana, A.A.; Hsu, A.; Beaver, L.M.; Wong, C.P.; Stevens, J.F.; Choi, J.; Jiang, Y.; Bella, D.; Williams, D.E.; et al. Untargeted Metabolomic Screen Reveals Changes in Human Plasma Metabolite Profiles Following Consumption of Fresh Broccoli Sprouts. Molecular Nutrition & Food Research 2018, 62, 1700665, doi:10.1002/mnfr.201700665.
  59. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B (Methodological) 1995, 57, 289–300.
  60. McMurdie, P.J.; Holmes, S. Phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. PLOS ONE 2013, 8, e61217, doi:10.1371/journal.pone.0061217.
  61. Oksanen, J.; Kindt, R.; Legendre, P.; Hara, B.; Simpson, G.; Solymos, P.; Henry, M.; Stevens, H.; Maintainer, H.; Oksanen@oulu,  jari The Vegan Package. 2009.
  62. Wickham, H. Ggplot2. WIREs Computational Statistics 2011, 3, 180–185, doi:10.1002/wics.147.
  63. Dunn, O.J. Multiple Comparisons Using Rank Sums. Technometrics 1964, 6, 241–252, doi:10.1080/00401706.1964.10490181.
  64. Dunn, O.J. Multiple Comparisons among Means. Journal of the American Statistical Association 1961, 56, 52–64, doi:10.1080/01621459.1961.10482090.
  65. Hartigan, J.A.; Wong, M.A. A K-Means Clustering Algorithm. Journal of the Royal Statistical Society: Series C (Applied Statistics) 1979, 28, 100–108, doi:10.2307/2346830.
  66. Segata, N.; Izard, J.; Waldron, L.; Gevers, D.; Miropolsky, L.; Garrett, W.S.; Huttenhower, C. Metagenomic Biomarker Discovery and Explanation. Genome Biology 2011, 12, R60, doi:10.1186/gb-2011-12-6-r60.
  67. Bates, D.; Mächler, M.; Bolker, B.; Walker, S. Fitting Linear Mixed-Effects Models Using Lme4. Journal of Statistical Software 2015, 67, 1–48, doi:10.18637/jss.v067.i01.
  68. Luang-In, V.; Albaser, A.A.; Rossiter, J.T. Characterization of a Recombinant β-Glucosidase of GH3 Family from Glucosinolate-Metabolizing Human Gut Bacterium Enterococcus Casseliflavus CP1 for Nitrile Production. Songklanakarin Journal of Science and Technology (SJST) 2020, 42, 549–556, doi:10.14456/sjst-psu.2020.69.
  69. Lai, R.-H.; Miller, M.J.; Jeffery, E. Glucoraphanin Hydrolysis by Microbiota in the Rat Cecum Results in Sulforaphane Absorption. Food Funct. 2010, 1, 161–166, doi:10.1039/C0FO00110D.
  70. Gil, V.; MacLeod, A.J. The Effects of PH on Glucosinolate Degradation by a Thioglucoside Glucohydrolase Preparation. Phytochemistry 1980, 19, 2547–2551, doi:10.1016/S0031-9422(00)83916-3.
  71. Agerbirk, N.; Olsen, C.E.; Sørensen, H. Initial and Final Products, Nitriles, and Ascorbigens Produced in Myrosinase-Catalyzed Hydrolysis of Indole Glucosinolates. J. Agric. Food Chem. 1998, 46, 1563–1571, doi:10.1021/jf9708498.
  72. Matusheski, N.V.; Swarup, R.; Juvik, J.A.; Mithen, R.; Bennett, M.; Jeffery, E.H. Epithiospecifier Protein from Broccoli (Brassica Oleracea L. Ssp. Italica) Inhibits Formation of the Anticancer Agent Sulforaphane. J. Agric. Food Chem. 2006, 54, 2069–2076, doi:10.1021/jf0525277.
  73. Williams, D.J.; Critchley, C.; Pun, S.; Nottingham, S.; O’Hare, T.J. Epithiospecifier Protein Activity in Broccoli: The Link between Terminal Alkenyl Glucosinolates and Sulphoraphane Nitrile. Phytochemistry 2008, 69, 2765–2773, doi:10.1016/j.phytochem.2008.09.018.
  74. Wathelet, J.-P.; Iori, R.; Leoni, O.; Rollin, P.; Mabon, N.; Marlier, M.; Palmieri, S. A Recombinant β-O-Glucosidase from Caldocellum Saccharolyticum to Hydrolyse Desulfo-Glucosinolates. Biotechnology Letters 2001, 23, 443–446, doi:10.1023/A:1010322322867.
  75. Lu, M.; Hashimoto, K.; Uda, Y. Rat Intestinal Microbiota Digest Desulfosinigrin to Form Allyl Cyanide and 1-Cyano-2,3-Epithiopropane. Food Research International 2011, 44, 1023–1028, doi:10.1016/j.foodres.2011.03.001.
  76. Luang-In, V.; Deeseenthum, S.; Udomwong, P.; Saengha, W.; Gregori, M. Formation of Sulforaphane and Iberin Products from Thai Cabbage Fermented by Myrosinase-Positive Bacteria. Molecules 2018, 23, 955, doi:10.3390/molecules23040955.
  77. Berteau, O.; Guillot, A.; Benjdia, A.; Rabot, S. A New Type of Bacterial Sulfatase Reveals a Novel Maturation Pathway in Prokaryotes*. Journal of Biological Chemistry 2006, 281, 22464–22470, doi:10.1074/jbc.M602504200.
  78. Charron, C.S.; Vinyard, B.T.; Ross, S.A.; Seifried, H.E.; Jeffery, E.H.; Novotny, J.A. Absorption and Metabolism of Isothiocyanates Formed from Broccoli Glucosinolates: Effects of BMI and Daily Consumption in a Randomised Clinical Trial. British Journal of Nutrition 2018, 120, 1370–1379, doi:10.1017/S0007114518002921.
  79. Wu, Y.; Shen, Y.; Zhu, Y.; Mupunga, J.; Zou, L.; Liu, C.; Liu, S.; Mao, J. Broccoli Ingestion Increases the Glucosinolate Hydrolysis Activity of Microbiota in the Mouse Gut. International Journal of Food Sciences and Nutrition 2019, 70, 585–594, doi:10.1080/09637486.2018.1554624.
  80. Bheemreddy, R.M.; Jeffery, E.H. The Metabolic Fate of Purified Glucoraphanin in F344 Rats. J. Agric. Food Chem. 2007, 55, 2861–2866, doi:10.1021/jf0633544.